def apply_bandpass_to_dataset(ds, variable_name='zos', lowcut=1/5, highcut=1/2, fs=24):
"""
Applies the band-pass filter to every pixel in an xarray Dataset.
"""
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
sos = butter(1, [low, high], btype='band', output='sos')
def _butter_wrapper(x):
return sosfiltfilt(sos, x)
filtered_ssh = xr.apply_ufunc(
_butter_wrapper,
ds[variable_name],
input_core_dims=[['time']],
output_core_dims=[['time']],
vectorize=True,
dask="parallelized",
output_dtypes=[ds[variable_name].dtype]
)
return filtered_ssh
ds_combined['ssh_filtered'] = apply_bandpass_to_dataset(ds_combined)Project A Update and Workflow
General Overview:
Hypothesis: Coastally trapped Internal Kelvin Waves (iKWs) are responsible for drastic swings in sea surface temperature (SST) through the vertical displacement of the thermocline - upwelling. We are trying to use modeled Sea Surface Height (SSH) data to characterise surface signals of iKWs around the coast of Newfoundland. We are then trying to identify corresponding signals of iKWs in water temperature data from moorings around the Newfoundland coastline.
GLORYS data:
For my project I am using GLORYS data, which is a global ocean reanalysis product. The variables I will be looking at are water temperature (WT) and sea surface height (SSH). The data is available on the Copernicus Marine Service website: https://marine.copernicus.eu/. This data is available at the temporal resolution of hourly for SSH and 6-hourly for WT, and at a spatial resolution of 1/12 degree. GLORYS is provided by the Mercator Ocean International consortium and is based on a combination of satellite observations, in-situ measurements, and ocean models. The data is available in NetCDF format, which can be easily read and processed using Python libraries such as xarray and netCDF4.
Station Data:
The other dataset I will be using in this analysis will be temperature data from the Headlands project. This data consists of 10 sites around the Newfoundland coastline where temperature loggers are deployed as moorings seasonally (see Figure 1). This data is available from 1989 to 2024. Then through NetCDF files, I can access water temperature data at the recorded depths and TimeStamp.
Water Temperature Analysis:
For the water temperature data, we construct a side-by-side time-series analysis of the mooring data and the GLORYS temperature data. This analysis gives insight into the performance of the GLORYS model in capturing the temperature variability at the mooring sites. We can also look for any correlations between the two datasets, which may indicate that the GLORYS model is successfully capturing the temperature swings (potential upwelling) signals associated with iKWs. Through this analysis we were able to move away from the HYCOM model, which was not performing well in capturing the temperature variability at the mooring sites, and instead focus all attention on the GLORYS model, which showed much better performance in this regard.
Based off of this analysis, we can then move forward with the GLORYS model and use it to identify potential iKW signals in the SSH data, and then look for corresponding signals in the temperature data at the mooring sites.
Sea Surface Height Analysis:
Identifying iKW signals in the SSH data is a bit more complex than the temperature data, as SSH is influenced by a variety of factors such as tides, barotropic waves, and ocean currents. To isolate the iKW signals, we can use a combination of filtering techniques and spectral analysis.
Temporal filtering:
We can apply a band-pass filter to the SSH data to isolate the frequency range associated with iKWs. This will help to remove any high-frequency noise and low-frequency signals that are not related to iKWs. Previous literature suggests that iKWs typically have periods ranging from 1 - 5 days.
Example code below using a Butterworth band-pass filter to isolate the iKW signals in the SSH data. The apply_bandpass_to_dataset function takes an xarray Dataset and applies the band-pass filter to the specified variable (default is ‘zos’ for sea surface height). The filter is designed to pass frequencies between lowcut and highcut, which are defined in terms of the sampling frequency fs. The filtered SSH data is then added to the original dataset as a new variable called ‘ssh_filtered’.
Spatial filtering:
For the spatial filter we can use a 2D rolling mean to calculate local SSH anomalies which will help to remove any large-scale signals that are not related to iKWs. This will help to isolate the local SSH signals that are more likely to be associated with iKWs. For this analysis we apply a window size of 12 x 12 which is approximately 1° of arc (or 100 km) in both latitude and longitude.
def calculate_local_anomaly(ds, variable='ssh_filtered', window_size=12):
"""
Calculates spatial anomaly using a NaN-aware rolling mean.
"""
local_mean = ds[variable].rolling(
latitude=window_size,
longitude=window_size,
center=True,
min_periods=1
).mean()
anomaly = ds[variable] - local_mean
ds_anomaly = ds.copy()
ds_anomaly[f'{variable}_anomaly'] = anomaly
return ds_anomaly
ds_processed = calculate_local_anomaly(ds=ds_combined, variable='ssh_filtered', window_size=12)It can be seen that applying the band-pass filter and spatial SSH anomaly filter separately does not fully isolate the iKW signals. Therefore, we can apply both filters to the SSH data to better isolate the iKW signals. This will help to remove any remaining noise and large-scale signals that are not related to iKWs.